Likelihood-free quantum inference: tomography without the Born Rule 
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We introduce the lilcelihood-free quantum inference algorithm (LFQIA), a Bayesian learning al- 
gorithm which can estimate quantum states and processes given access to a classical experiment 
simulator rather than to exact likelihood functions. It has been shown recently that there exists 
interesting classes of states and processes for which weak simulation (the ability to efficiently sample 
the probability distribution of outcomes of the experiment) is possible but strong simulation (the 
ability to efficiently compute the full distribution of outcomes of the experiement) is not. In such 
cases, our algorithm exponentially outperforms standard algorithms based on computation of the 
likelihood function and we demonstrate this with an example. 



Practical quantum information processing requires 
precise control over physical systems that store and ma- 
nipulate quantum information. Characterizing quantum 
information processing devices is required for develop- 
ing, testing, tuning and, ultimately, controlling them. 
This characterization process is generically called quan- 
tum state tomography or quantum, process tomography [1] . 
Though tomographic experiments have been carried out 
in systems as large as several qubits [2, ^ , the number of 
parameters that must be measured in tomography grows 
exponentially with the number of qubits. As a result, 
in order to practically characterize systems large enough 
to be useful in quantum information processing appli- 
cations, devising novel and efficient statistical inference 
methods is a necessity. 

There already exist many ingenious methods aimed at 
reducing the complexity of identifying and characteriz- 
ing quantum devices. These include identifying stabi- 
lizer states 0] ; tomography for matrix product states [5] ; 
learning local Hamiltonians |^ ; tomography for low-rank 
states via compressed sensing [7]; and tomography for 
multi-scale entangled states [8] . These technqiues employ 
efficient simulation algorithms which propagate a repre- 
sentation of the state vector and provide a calculation of 
the probabilities of the distribution of measurement out- 
comes given by the Born rule. When such a calculation is 
efficient, it is coined strong simulation [9Ullj. Thus, pro- 
vided the model under study is such that classical strong 
simulation can be performed efficiently, estimation can 
also be done efficiently. However, if the state or process 
in question is not efficiently simulatable, the estimation 
algorithm will also not be efficient. 

More recently there has been growing interest in 
an alternative and less restrictive notion of simulation. 
Namely, weak simulation requires only that the algo- 
rithm efficiently produce samples from the output prob- 
ability distribution [9l[TT|. Validation techniques (esti- 
mating "closeness" to some target) have recently been 
proposed which take advantage of this kind of simula- 



tion [6l [T2l415j . In this paper, we provide a method for 
state and process estimation when strong simulation is 
not possible and without explicit calculation of the prob- 
ability distribution of outcomes (the likelihood function) . 
In such cases, our algorithm exponentially outperforms 
standard statistical techniques. 

Since our algorithm needs only a classical simulator 
and does not require calculation of the likelihood func- 
tion itself, we call our algorithm the likelihood-free quan- 
tum inference algorithm (LFQIA). For brevity, we will 
phrase the general problem of statistical inference in the 
language of quantum state estimation. However, it is cru- 
cially important to remember that everything said will be 
true for estimating quantum processes, measurements, or 
any other quantity that can be parameterized by a set of 
real numbers and for which efficient simulation, in the 
weak sense, is available. 

For state estimation, then, the task is as follows. Given 
a set of data specified as a set of observable-eigenvalue 
pairs D := {Ok,ek} and assumed to be measured on 
copies of the same state UpU^ for some known evolution 
operator U, we are to estimate the unknown initial state 
p. All probabilistic methods feature the likelihood func- 
tion, the Born rule, which specifies the probability of this 
data set, given some state UpW: 



PTiD\p;U)=Y[{Ok\UpU^\Ok) 



(1) 



where \Ok) is the eigenvector of ek in the basis diago- 
nalizing Ok- The current standard is to maximize this 
function, selecting as our estimated state the one which 
achieve this maximum: the maximum likelihood estima- 
tor [TB]. Some would refer to this as, in the language of 
interpretations of probability, frequentism. 

On the other end of the spectrum is Bayesianism. In 
this approach, one first chooses a prior probability dis- 
tribution Pr(yo) over the set quantum states. Then, via 
Bayes' rule, we arrive at the distribution over states given 



the observed data (the posterior distribution) 



Pt{p\D:U) = 



Pr{D\p;U)PT{p) 
JPT{D\p;U)PT{p)dp- 



(2) 



For a staunch Bayesian, we are done; the posterior rep- 
resents complete knowledge of the parameters. We, how- 
ever, will demand the more practical and efficient require- 
ment of specifying an estimate of the mean and covari- 
ance of the posterior distribution of unknown parameters. 
Bayesian estimation of this kind was first suggested for 
application to quantum tomography by Blume-Kohout 

Having adopted the Bayesian approach, then, the con- 
crete problem is to provide an efficient numerical algo- 
rithm performing Bayesian inference in the above sense. 
Moreover, we require this without a means to calculate 
the likelihood function (requiring only a means to sample 
from it). Next, we give an overview of the algorithm we 
propose to solve this problem. 

In order to efficiently compute the integral in equa- 
tion (I2]), we employ the sequential Monte Carlo (SMC) 
method, which has been used for the purpose of Hamil- 
tonian learning [18' and in the tomographic estimation of 
one and two qubit states |19j , as well as in the continuous 
measurement of a qubit [SO] . 

The sequential Monte Carlo method prescribes that 
we should approximate a distribution over model param- 
eters with a distribution that has support only over a 
finite number of points (often referred to as particles). 
Each particle is assigned a weight (colloquially thought 
of as it's relative plausibility). More concretely, we ap- 
proximate an arbitrary distribution by 



Pr(p|Z?)«^Wfe(i?)%-pfc) 



(3) 



fc=i 



where the weights at each step are iteratively calculated 
from the previous step via 



Wkid-j+i) oc PT{dj+i\pk)wk{dj), 



(4) 



followed by normalization. Using this form of a sequen- 
tial Monte Carlo algorithm also has the advantage of en- 
suring that the prior and posterior distributions for the 
update always have support over a finite number of par- 
ticles, which simplifies the analysis of our algorithm. We 
explored some variants of this algorithm and presented 
it in much greater detail in reference |18| . 

The particle approximation can be made arbitrarily 
accurate by increasing the number of particles and will 
be a good approximation at every update provided we 
feed in, at the initial stage, the appropriate weights {wk} 
and support points {pk}- Such a choice is to have equal 
weight, Wk — 1/n for all k, and sampled particles accord- 
ing to the correct prior Pr(p). 

So far, we seem to require a full specification of the like- 
lihood function Pr(dfe|p). To free us from this crutch, we 



reconstruct the likelihood function from simulated data. 
Recall, from the definition, when we have access to a 
weak simulator, the step of sampling the likelihood func- 
tion will be efficient. For brevity, we will discuss the bi- 
nary case with outcomes labeled and 1. The unknown 
probability po := Pr(0|p) can be treated as a parameter 
to be estimated. In particular, since we have assumed 
that the data are conditionally independent given the 
model, repeatedly sampling the likelihood function will 
produce data that follows a binomial distribution with 
parameter pq. Estimating the parameter of a binomial 
distribution from sample data is a well-understood sta- 
tistical problem. Supposing k O's were observed in n 
trials, a typical estimator is 



POf{k) = 



k + 'j 

n + 2j' 



(5) 



where 7 is a free parameter. These are called "linear" or 
"add-7" estimators [H]. The latter phrase is due to the 
equivalence to standard maximum likelihood estimation 
when adding 7 "fake" observations — also termed "hedg- 
ing" [22]. These estimators can also be understood to 
arise from a Bayesian approach as well. In particular, 
the estimator in equation ([5]) is the posterior mean when 
using the following Beta distribution as a prior 



Pr(po) ocPo ^{l~pa) 



7-1 



(6) 



The posterior variance of this distribution can also be 
calculated as 



3.2 (i^-. ^ {k + -/){n - k + j) ^ Po^jl-Poj) 
^°'^^ ' (n + 27)2(n + 27 + l) n -t- 27 + 1 ' 



(7) 



Here we will use the value 7 = 1, as it corresponds to a 
uniform prior distribution. We leave the optimization of 
this algorithmic parameter for future work. 

If we are willing to tolerate an error e in our recon- 
struction of the likelihood, then we can check after each 
sample if ctq^ < e. If not, we collect more samples until 
the condition is met. We therefore have a single quality 
parameter for this adaptive protocol: e. 

We illustrate the utility of LFQIA by applying it to 
an example tomographic problem for which it is hard 
to compute the likelihood function but easy to simulate 
from it. Introduced by Van den Nest [10], the example 
involves a class of circuits which can be efficiently sim- 
ulated in the weak sense but is hard to simulate in the 
strong sense. The family of circuits are those that ap- 
ply Hadamard gates followed by "classical" gates (NOT, 
CNOT and Toffoli) and measure the first qubit in the 
computational basis. These are called HT-circuits, and 
are depicted in Figure [T] 

The input state is the "all-zero" state |0) . If a 
Hadamard gate is applied to each of the first m of 
the n total qubits, and if the classical gates compute 



|0) -{H} 

|o) ^0- 
|o) ^0- 

{Pk,Pi,} — 



V^ 



Uf 



where Fy{x) = f{x,y). Note the original situation is the 
special case of Fq. The likelihood function (the probabil- 
ity of data given the distribution py) is now 



Pr(d|/,p,) = ^p, 



\{xe{0,ir':F,{x)i^d}\ 



(10) 



FIG. 1: An HT circuit. There are m qubits which are acted 
upon by Hadamard gates and n qubits in total. In their 
original construction, demonstrating the distinction between 
strong and weak simulation [TU], the final n — m qubits also 
are initialized in the state ]0). This becomes a tomography 
problem if we allow these ancilla states to be unknown. 



a polynomial-time computable function / : {0, 1}" — > 
{0, 1}", then the state evolves to 

|0)®" ^ ;7/if^" |0)®" = Y, I^W), (8) 

where F{x) := f{x, 0, . . . , 0). The measurement yields a 
bit d E {0, 1} with probability 

Pr(d|F)='^-^^»'^>;;^^-^^--^>', (9) 

and now F{x)i denotes the first bit of F{x). 

To generate a sample from this distribution is easy: 
first sample an m-bit string x uniformly at random, then 
compute F{x) and set d = F{x)i. On the other hand, 
computing this probability is hard since, in order to do 
so for an arbitrary F, we must evaluate F on all inputs. 
Given that there are 2™ inputs, this problem is exponen- 
tially hard. More precisely, the problem is ttP-complete, 
as can be seen by considering a polynomial-time Turing 
machine M and letting Jm '■ {0, 1}" -^ {0, 1}" repre- 
sent M by assigning fM{x) = 0" if M accepts x and 
fnix) = 1" otherwise. Evaluating the likelihood func- 
tion for the circuit when / — /m is thus equivalent to 
counting the accepting paths of M jJTj. Thus, our de- 
mand for arbitrary precision in the calculation of the like- 
lihood function has caused us to have to consider all pos- 
sible paths of a polynomial-time Turing machine, while 
sampling demands only that we sample one path per sam- 
ple. 

This becomes a tomography problem if we suppose the 
last n—m qubits (the ones not acted upon by Hadamards) 
are in unknown computational basis states (see again Fig- 
ure fl]) . More precisely, suppose the state that is input 
is drawn from the unknown ensemble {py, \y)}, where 
y G {0,1}"^™. The goal is to estimate the distribution 

Py 

In this more general setting, given y, the final state is 



Since this situation is more general than the special case 
considered above, the same argument for the hardness 
of strong simulation applies. It is still easy, however, to 
simulate outcomes from this distribution. Indeed, the 
exact same method as before works, with one small addi- 
tion: first, draw y from py] then draw an m-bit string 
X uniformly at random; then compute Fy{x) and set 
d = Fy{x)i. 

We specialize this example to the case where the ancilla 
state is described by the ensemble density matrix p = 
p|0)(o|»"-™ -f (1 - p)|l)(lp"-™. This example is ideal 
for illustrative and testing purposes because it contains 
enough complexity to make the technique relevant but 
can also be solved analytically. 

In estimating the probability p, we benchmark LFQIA 
using the mean squared error E = Ep[Efe|p[|p — p(fc)| ], 
where p{k) is the estimator produced by our algorithm 
after observing k Os at the output of the HT circuit. We 
compare LFQIA to our previous SMC algorithm (which 
calculates the likelihood function exactly), the exact so- 
lution and the following asymptotic expression we have 
derived for the mean squared error: 



E = 



1 



6(po-Pi)' 



1 

TV 



1 

iV2 



o 



1 

iV3 



(11) 



xG{0,l} 



Fy[x)), 



where po = Pi'(0|/, |00 . . . 0)) and where pi is defined sim- 
ilarly. Both po and pi can be computed offline. 

The performance of LFQIA for this problem is illus- 
trated in Figure l2] and, taking computational costs into 
account. Figure |3j Figure [2] shows the error in estima- 
tion for the exact solution as well as for LFQIA. As 
expected, LFQIA — having additional approximations — 
performs less favorably than the rest. However, there 
is strong evidence that the algorithm is consistent in 
the sense that, as the tolerance decreases, the perfor- 
mance increases. This conclusion has been reached inde- 
pendently by verifying the posterior distributions at the 
output are close in a metric on the space of distributions 
(these results will be presented in detail elsewhere). 

As shown in Figure |3J when one considers the com- 
putational cost of the algorithm, the conclusion above 
is reversed: LFQIA is now optimal when the number 
of qubits renders computing the likelihood function in- 
feasible. The performance of LFQIA illustrates a clear 
trade-off in accuracy versus efficient. However, the rela- 
tively small penalty in accuracy is overwhelmed by the 
exponential improvement in efficiency as the system size 



mcreases. 
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FIG. 2: The mean squared error of the estimator produced 
by the exact solution as well as SMC and LFQIA. The prior 
for this example was chosen to be uniform and thus the ex- 
act solution (red line) was tractable. The asymptotic mean 
squared error (black line) is independent of the prior and given 
by equation ( |11[ ). The remaining two lines are the actual per- 
formance of SMC and LFQIA. In all cases, n = 100 particles 
were used. As expected, the exact solution is best in all cases. 
The SMC method should also perform better than LFQIA 
when the same number of particles are used, as we see. As 
the error tolerance parameter of LFQIA, e, decreases, the ac- 
curacy improves. Since this comes at the cost of requiring 
more samples, we exhibit a clear trade-off in efficiently and 
accuracy. 



In this work, we have demonstrated an improvement 
of the sequential Monte Carlo parameter estimation al- 
gorithm that allows for its extension to the case of weak 
(sampling) simulators. For models with fast weak sim- 
ulation available, our algorithm can be seen to provide 
dramatic advantages in terms of classical computing costs 
over sequential Monte Carlo alone and at minimal cost 
in estimation performance. This extension allows for us 
to perform inference in subtheories of quantum mechan- 
ics that admit a large separation between the tractibility 
of strong and weak simulation. In this way, our algo- 
rithm demonstrates more concretely the deep connection 
between simulating a model and estimating properties of 
that model. As the complexity of candidate quantum in- 
formation processors grows, our algorithm provides a way 
forward to estimating properties of very large systems by 
exploiting this deep connection. Moreover, our results 
demonstrate that LFQIA provides an approximately cor- 
rect representation of the posterior distribution such that 
quality parameters may be chosen according to practical 
needs and available resources. 

We have focused here on the HT-circuit model for its 
conceptual simplicity and analytical convenience. How- 
ever, we expect similar findings for other scenarios where 
weak simulation is used. In particular, a large class of 
such schemes have been proposed for simulating quan- 
tum dynamics which constrain the Wigner function of 
the quantum state to be positive [23H25] . Moreover, 
these method could also be generalized to other quasi- 
probability distributions — all of which demonstrate that 
the range of applicability and the enhancement provided 
by LFQIA could be quite large indeed. 
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FIG. 3: The average reduction in mean squared error of 
the estimator per second of computation time produced by 
SMC and LFQIA. The prior for this example was chosen to 
be uniform. In all cases, n — 100 particles were used and the 
computation time and mean squared error was recorded after 
N = 100 measurements. What the figures are meant to high- 
light and clearly show is the trade-off in accuracy and com- 
putation time for both SMC and LFQIA. The most impor- 
tant message, however, is that SMC, requiring evaluations of 
the likelihood, is exponentially hard in the number of qubits. 
Thus, per unit of computation time, SMC is horribly ineffec- 
tive as an estimation algorithm for such problems. On the 
other hand, LFQIA achieves the same level of performance 
independent of the number of qubits. The reason that higher 
error tolerance, e, seems to perform better in this parameter 
regime is due to the fact that higher tolerance requires fewer 
samples and hence less computation time. Again, this demon- 
strates a trade-off in the acciuracy versus efficient in learning 
parameters in quantum mechanical models. 
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